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ABSTRACT 


We  are  continuing  efforts  to  develop  3-D  velocity  models  for  southern  Asia  through  the  collection  and  analysis 
of  broadband  waveform  data  acquired  on  the  Indian  subcontinent.  The  geology  of  India  is  diverse,  but  can  be 
divided  into  three  main  regions:  the  Himalayas,  the  Indo-Gangetic  plain,  and  the  Indian  Shield.  Our  initial 
focus  has  been  on  the  southern  Indian  shield,  and  we  have  also  completed  studies  in  the  trans-Himalayas  and 
the  Shillong  Plateau.  The  goal  of  our  work  is  to  determine  the  crust  and  upper  mantle  velocity  and  attenuation 
structure  and  to  characterize  regional  seismic  waveform  propagation  of  the  Indian  subcontinent.  Teleseismic 
receiver  function  data,  S-to-P  conversions,  and  short-period  surface  wave  phase  velocity  data  have  been 
interpreted  for  seismograms  recorded  along  a  700-km  north-south  profile  of  the  southern  Shield.  These  data 
show  that  the  shield  velocity  structure  is  extremely  uniform,  simple,  and  consists  of  a  surface  wave  velocity  of 
approximately  3.45  km/sec  and  a  moderate  gradient  of  0.20  km/sec/km  with  the  Moho  at  35±1  km  depth.  To 
the  south  of  the  shield  in  the  granulite  terrain,  the  crust  is  both  thicker  (44+ 1  km  deep)  and  more  complicated, 
with  a  mid-crustal  discontinuity  at  approximately  25-km  depth.  We  constrain  the  upper  mantle  structure  with 
phase  velocity  measurements  of  long-period  surface  waves;  these  show  that  the  seismic  lithosphere  is 
approximately  150  km  thick  and  underlain  by  a  weak  low-velocity  zone.  We  have  also  modeled  the  regional 
waveforms  of  a  moderate  earthquake  that  occurred  near  Koyna,  India,  in  September  2000  and  of  the  26  January 
200 1  Bhuj  main  shock  and  its  aftershocks  to  calibrate  paths  to  the  regional  stations. 
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OBJECTIVES 


The  objectives  of  our  research  effort  are  fourfold.  The  first  is  to  increase  the  amount  of  high  quality  seismic  data 
that  are  available  by  deploying  additional  seismic  stations  in  the  India/Pakistan  region.  The  second  is  to  use 
these  data  to  refine  the  regional  attenuation  characteristics,  the  regional  travel  time  corrections,  and  the  crust 
and  upper  mantle  structure  of  the  region.  The  third  is  to  use  the  analysis  of  these  data  to  refine  Weston's  3-D 
velocity  model  for  India  and  Pakistan  (WINPAK3D)  and  to  validate  the  model.  The  fourth  is  to  provide  a 
cooperative  forum  for  the  exchange  of  information  and  ideas  on  the  analysis  of  the  data  and  modeling  of  the 
region. 

This  research  is  motivated  by  the  fact  that  calibrated  regional  velocity  models  for  the  India/Pakistan  region  will 
be  essential  to  improve  upon  the  location  capability  of  global  networks.  For  example.  Table  1  lists  five 
earthquakes  that  occurred  in  the  Koyna-Warna  region  along  the  western  coast  of  India.  Yet,  even  with  some 
overlap  in  the  datasets  processed  by  the  three  institutions,  the  locations  formed  by  the  networks  show  significant 
scatter  (Figure  1).  In  fact,  of  the  15  different  locations  for  the  five  events,  only  four  locations  fall  within  or  near 
the  boundary  of  the  seismicity  cluster  regions  determined  for  the  Koyna-Warna  region  by  Mandal  et  al.  (1998). 
This  region  is  well  known  for  the  clustered  seismic  activity  that  started  in  1962  after  the  impoundment  of  the 
Koyna  reservoir.  The  events  were  recorded  on  global  seismic  networks  and  located  by  the  Prototype 
International  Data  Center  (pIDC),  United  States  Geological  Survey  (USGS),  and  International  Seismic  Centre 
(ISC).  In  addition  to  the  locations,  we  present  the  range  of  group  velocity  dispersion  curves  generated  for  HYB 
data  (Figure  2)  using  the  event  source  locations  and  origin  times  presented  in  Table  1.  For  reference,  we  have 
plotted  the  theoretical  dispersion  curve  generated  using  our  shear  wave  velocity  model  for  the  southern  Indian 
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shield  that  is  discussed  later  in  this  paper  and  presented  in  Figure  4.  Because  the  path  from  these  events  to  HYB 
(A  -  525  km;  azimuth  -  90°)  is  essentially  the  same  for  all  events  in  the  region,  and  focal  mechanisms  for  all 
three  clusters  are  similar  (Mandal,  et  al,  1998),  differences  in  the  dispersion  curves  must  be  related  to  location 
and/or  origin  time  errors.  The  average  dispersion  curve  generated  by  these  data  at  shorter  periods  has 
approximately  ±  0.2  seconds  error  and  thus  would  incorporate  considerable  error  if  used  to  develop  a  shear 
wave  velocity  model  for  the  region.  For  the  model  development  described  in  the  following  sections  of  this 
paper,  we  have  searched  for  and  examined  only  well-located  events  in  order  to  reduce  such  model  error.  The 
results  are  high-quality  regional  velocity  models  for  southern  Asia  that  have  been  incorporated  into 
WINPAK3D  and  are  being  used  to  more  accurately  locate  regional  seismic  events  in  southern  Asia. 
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Table  I.  Source  parameters  from  the  pIDC,  ISC,  and  USGS  for  five  earthquakes  near  the  Koyna-Warna  region 
of  western  India. 
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Figure  1.  USGS,  pIDC,  and  ISC  network  locations  for  the  five  earthquakes  listed  in  Table  1.  The  three 
ellipsoids  represent  the  distinct  regions  of  the  Koyna-Warna  seismic  zone  where  the  majority  of  the 
seismicity  occurs  (Mandal  et  al.,  1998)  The  abbreviations  for  each  cluster  include  NK  -  Northern 
Koyna;  SK  -  Southern  Koyna;  and  B  -  Bhogiv. 
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Figure  2.  Group  velocity  dispersion  for  paths  between  the  Koyna-Warna  region  and  station  HYB.  The  scatter 
in  the  dispersion  curves  is  related  to  location  and  origin  time  errors  and  illustrates  the  need  for  using  only 
high-quality  event  data  to  generate  the  calibrated  regional  velocity  models.  The  blue  line  is  the 
theoretical  dispersion  curve  for  the  south  Indian  shield  as  determined  by  research  discussed  in  the 
following  section. 


RESEARCH  ACCOMPLISHED 


Installation  of  Broadband  Seismographs  in  India 

In  early  February  we  received  five  broadband  seismographs  from  Guralp  Systems  Limited,  which  we  deployed 
in  India  in  March/ April  of  this  year.  Figure  3  shows  the  location  of  these  stations  and  other  stations  we  are 
currently  operating  in  India  in  cooperation  with  the  National  Geophysical  Research  Institute  at  Hyderabad  and 
the  Indian  Institute  of  Astrophysics  at  Bangalore. 

Crustal  structure  beneath  southern  India 


Short-period  surface  wave  studies.  Surface  wave  phase  velocities  for  fundamental  mode  Rayleigh  waves 
have  been  measured  for  a  number  of  paths  across  the  south  Indian  Shield  (Figure  3).  The  analysis  of  surface 
wave  dispersion  data  gives  an  average  crustal  model  for  the  south  Indian  Shield  that  is  used  as  a  starting  model 
in  the  receiver  function  inversions  discussed  below,  and  provides  additional  constraints  in  the  joint  receiver 
function  -  surface  wave  inversion. 

We  determined  an  average  short  period  fundamental  mode  Rayleigh  wave  dispersion  curve  using  various  two- 
station  pairs  along  the  North-South  (NS)  profile  using  the  transfer  function  method  of  Gomberg  et  al.  (1988). 
This  technique  poses  the  problem  of  phase  velocity  determination  as  a  linear  filter  estimation  problem  in  which 
the  seismogram  at  a  more  distant  station  from  an  event  is  the  convolution  of  the  seismogram  at  a  close  station 
with  the  Earth  filter  (the  phase  velocity  curve)  to  be  determined.  Smoothness  constraints  are  imposed  based  on 
an  approximate  knowledge  of  the  group  velocity.  We  used  the  dispersion  results  of  Bhattacharya  (1992)  as  an 
initial  dispersion  model.  To  ensure  the  starting  dispersion  curve  was  appropriate,  we  tested  a  number  of  initial 
dispersion  models  and  smoothing  criteria;  the  final  dispersion  curves  (Figure  4a)  are  the  ones  where  the  results 
were  stable  in  the  sense  that  the  resulting  dispersion  curve 
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Figure  3.  Map  showing  the  stations  locations:  triangles  -  permanent  stations  of  Indian  Institute  of 
Astrophysics/University  of  Cambridge;  dots  -  temporary  recording  sites  of  National  Geophysical 
Research  Institute/University  of  Cambridge;  squares  -  global  digital  seismic  network  stations. 


Shear  Wave  Velocity  (km/s) 


Figure  4.  Right  -  Fundamental  mode  Rayleigh  wave  dispersion  curves:  solid  line  with  error  bars  -  average 

dispersion  for  the  NS  profile;  dashed  lines  -  dispersion  bounds  for  the  northern  portion  of  NS  profile;  dot- 
dash  lines  -  dispersion  bounds  for  the  southern  portion  of  NS  profile.  Left  -  Crustal  velocity  model  derived 
from  inversion  of  the  average  dispersion  curve. 
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was  not  strongly  influenced  by  realistic  changes  in  the  initial  dispersion  model  and  smoothing  criteria.  The 
dispersion  curves  in  Figure  4a  were  determined  simultaneously  from  multiple  station  pairs,  but  we  also 
determined  dispersion  curves  for  all  individual  two-station  paths  separately  to  verify  that  there  were  no  large 
outliers  amongst  the  various  two-station  combinations. 

We  invert  the  Rayleigh  wave  phase  velocity  data  using  the  stochastic  least-squares  routine  of  Hermann  (1994). 
This  expresses  the  standard  least-squares  problem  in  terms  of  eigenvalues  and  eigenvectors  and  uses  singular 
value  decomposition  to  invert  the  matrix  giving  the  solution  vector,  the  variance-covariance  matrix,  and  the 
resolution  matrix.  Since  the  dispersion  data  for  the  northern  portion,  southern  portion,  and  whole  profile  were 
nearly  identical,  and  since  the  crustal  models  from  refraction  data  at  the  northern  (Krishna  et  al.,  1999)  and 
southern  (Kr  ishna  and  Ramesh,  2000)  ends  of  the  profile  are  similar,  we  inverted  only  the  average  dispersion 
curve.  The  starting  model  for  the  inversion  is  from  the  average  crustal  model  for  the  Indian  Shield  of  Singh  et 
al.  (1999)  and  consisted  of  an  upper  crustal  layer  13.8  km  thick  with  Vs  3.55  km/s  and  a  lower  layer  24.9  km 
thick  with  Vs  3.85  km/s,  over  a  4.65  km/s  half-space.  The  final  inversion  model  and  the  fit  of  the  dispersion 
curve  for  this  model  to  the  observed  dispersion  are  shown  in  Figure  4.  The  average  crustal  structure  beneath  the 
profile  consists  of  two  layers  (hi  12  km,  Vsj  3.65  km/s,  h2  23  km,  Vs2  3.81  km/s)  overlying  a  mantle  with  S„ 
velocity  4.61  km/s. 

Receiver  function  analysis.  Figures  5  and  6  illustrate  the  receiver  function  analysis  we  have  performed  using 
data  from  Bangalore  (BGL).  BGL  is  situated  on  the  western  edge  of  the  eastern  Dharwar  craton  close  to  its 
contact  with  the  Closepet  granitic  intrusion  (Figure  3).  Figure  5  shows  the  BGL  receiver  function  with  respect 
to  back  azimuth.  The  eastern  quadrant  is  well  sampled,  but  there  is  only  one  event  from  the  west.  The  small 
variation  in  the  Ps  converted  phase  (4.0±0.1)  and  the  small  amplitude  of  the  tangential  receiver  function  relative 
to  the  amplitude  of  the  radial  receiver  function  justify  modeling  the  crust  beneath  BGL  as  1-D. 

Figure  6  (left)  shows  a  receiver  function  stack  for  BGL  consisting  of  10  events  with  distance  76±7°  and  back 
azimuth  99±3°  E.  The  receiver  function  has  a  large  amplitude  positive  arrival  at  about  4  sec  (PS)  followed  by  a 
prominent  positive  arrival  at  about  14  seconds  (PpPmS)  and  a  negative  arrival  at  16-18  seconds  (PpSmS  + 
PSPmS)).  These  phases  are  also  apparent  in  many  of  the  individual  receiver  functions  making  up  the  stack. 
Below  the  stack  is  the  tangential  receiver  function.  Amplitudes  of  the  tangential  arrivals  are  small  compared  to 
those  of  the  radial  arrivals. 

The  receiver  function  stack  was  inverted  for  1-D  crustal  velocity  models.  To  test  the  sensitivity  of  the  inversion 
to  the  starting  model,  we  used  three  different  initial  models:  (1)  the  refraction  model  for  the  GBA  region  of 
Krishna  and  Ramesh  (2000),  (2)  the  average  crustal  model  derived  from  the  surface  wave  phase  velocity 
dispersion  (Figure  4),  and  (3)  a  model  derived  from  the  timing  of  the  conversions  and  reverberations  in  the 
receiver  functions  and  the  average  P  wave  velocity  of  the  crust  from  Krishna  and  Ramesh  (2000).  Zandt  et  al. 
(1995)  pointed  out  that  the  average  crustal  properties  can  be  estimated  directly  from  the  arrival  time  differences 
of  the  pP,  PS,  PpPmS  and  PpSmS  +  PSPmS  phases  in  the  receiver  function.  The  PS  -  pP  time  difference  is 
dependent  on  the  average  Vp/Vs  ratio  of  the  crust  and  the  crustal  thickness.  The  PS  -  PpPmS  time  is  the  two-way 
P  wave  travel  time  and  the  {PpSmS  +  PSPraS}  -  pP  time  is  the  two-way  S  wave  travel  time  through  the  crust  for 
a  ray  with  ray  parameter  p.  The  ratio  of  PS  -  pP  to  PS  -  PpPmS  is  independent  of  crustal  thickness  but  weakly 
dependent  on  Vp;  the  ratio  of  PS  -  PpPmS  to  {PpSmS  +  PSPmS}  is  proportional  to  the  Vp/Vs  ratio  and 
independent  of  crustal  thickness. 

We  inverted  the  receiver  functions  using  the  linearized  inversion  of  Ammon  et  al.  (1990)  both  with  and  without 
the  surface  wave  dispersion  constraint  and  the  final  inversion  results  were  nearly  identical  for  the  three  starting 
models.  We  then  simplified  the  model  by  grouping  adjacent  model  layers  with  similar  wave  speeds  to  form  a 
more  coarsely  parameterized  starting  model  and  re  inverted  the  receiver  function.  We  repeated  this  procedure 
until  we  found  the  velocity  model  with  the  minimum  number  of  parameters  that  fit  both  the  main  features  of  the 
receiver  function  and  the  short  period  surface  wave  phase  velocity.  Finally,  we  tested  the  main  features  of  the 
crustal  model  (e.g.,  thickness  of  the  low  velocity  surface  layer,  thickness  of  the  gradient  at  the  base  of  the  crust, 
Moho  depth)  using  a  forward-modeling  model  to  estimate  how  well  these  model  features  were  constrained  by 
arrivals  in  the  observed  receiver  function. 
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BGL  ozimuthol  plot 


Figure  5.  BGL  receiver  functions  vs.  azimuth:  heavy  solid  line  -  radial  receiver  function;  red  solid  line  - 

tangential  receiver  function;  dashed  line  -  radial  receiver  function  calculated  for  the  final  crustal  velocity 
model  shown  in  Figure  6. 
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Figure  6.  Receiver  function  analysis  summary  for  BGL.  Right:  (a)  initial  inversion  results;  (b)  results  of  joint 
receiver  function/surface  wave  inversion;  (c)  final  inversion  models;  (d)  forward  modeling  test;  (e) 
match  between  the  observed  phase  velocity  and  the  phase  velocity  calculated  from  the  models  shown  in 
(b)  and  (c).  Left:  Match  of  the  receiver  function  calculated  for  the  four  velocity  models  at  right  and  the 
+  1  ct  bounds  for  the  BGL  receiver  function  stack. 
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These  steps  are  summarized  in  Figure  6  (right).  The  initial  inversion  model  without  the  dispersion  constraint 
shows  a  smooth  velocity  gradient  from  about  Vs  3.1  km/s  at  the  surface  to  about  3.75  km/s  in  the  lower  crust 
and  a  Moho  discontinuity  at  35±2  km  depth  (Figure  6a).  The  inversion  with  the  surface  wave  dispersion 
constraint  (Figure  6b)  produces  a  model  identical  to  those  shown  in  Figure  6a.  The  final  velocity  model  (Figure 
6c)  shows  a  thin  (2-4  km)  low  velocity  surface  layer,  an  almost  uniform  velocity  crust  (3. 5-3. 6  km/s),  and  a 
Moho  discontinuity  at  35±1  km  depth.  The  fits  of  the  synthesized  receiver  functions  for  the  various  velocity 
models  to  the  ±  1  a  bounds  of  the  observed  receiver  function  stack  are  shown  in  Figure  6  (left). 

In  the  forward  modeling  test,  the  upper  mantle  Pn  and  Sn  velocity  were  fixed  at  8.2  km/s  (Krishna  and  Ramesh, 
2000)  and  4.72  km/s  (Fluestis  el  al.,  1973),  respectively,  but,  in  fact  the  receiver  functions  are  relatively 
insensitive  to  these  parameters.  Neither  the  difference  in  the  near  surface  layer,  the  mid-crustal  discontinuity, 
nor  the  gradient  at  the  base  of  the  crust  from  the  inversion  of  the  different  stacked  receiver  functions  is 
significant,  and  we  find  that  considering  the  noise  in  the  data,  one  crustal  velocity  model  fits  all  the  receiver 
function  stacks  for  BGL.  Figure  5  compares  the  synthetic  receiver  function  for  the  BGL  crustal  model  which  the 
observed  receiver  functions  as  a  function  of  azimuth.  This  agreement  indicates  little  lateral  variation  in  the 
crustal  velocity  structure  about  BGL.  If  the  small  systematic  variation  in  the  PS  -  P  time  suggested  by  Figure  5 
is  real,  this  could  result  from  a  systematic  variation  of  crustal  velocity  about  BGL  or  a  systematic  variation  in 
Moho  depth.  If  the  crustal  velocity  variation  is  spread  through  the  whole  35  km  thickness  of  the  crust,  it  would 
correspond  to  a  ±0.1  km/s  variation  in  shear  wave  velocity.  If  the  variation  is  due  to  a  systematic  azimuthal 
variation  in  the  Moho  depth,  this  corresponds  to  a  ±1  km  lateral  variation  in  the  Moho  depth  about  BGL. 

We  have  analyzed  data  from  14  broadband  sites  in  southern  India  in  this  manner  and  the  resulting  Moho  depths 
are  shown  in  Figure  7.  The  Moho  is  at  almost  constant  depth  along  the  western  edge  of  the  eastern  Dharwar 
craton.  However,  the  two  sites  on  the  western  Dharwar  craton  show  the  Moho  there  is  about  5  km  deeper,  a 
result  also  noted  in  the  Kavali  -  Uduppi  refraction  study  (Kaila  and  Krishna.  1992).  To  the  south  of  the  shield  in 
the  granulite  terrain  (KOD,  Fig.  1)  the  crust  is  43  km  thick.  Thus,  there  are  significant  differences  in  Moho 
depth  over  southern  India,  but  the  results  we  have  to  date  suggest  they  are  well  correlated  with  the  surface 
geology. 

Upper  mantle  structure  beneath  southern  India.  Figure  8  shows  the  fundamental  mode  Rayleigh  waves 
dispersion  measurement  from  seven  two-station  pairs  using  the  combination  of  the  stations  NND/HYB  in  the 
north  and  GBA/KOD  in  the  south.  Inter-station  path  lengths  for  the  two-station  pairs  vary  from  600  to  800  km. 
The  coherency  plot  shows  that  the  dispersion  is  well  constrained  to  — 180  sec  period.  Previous  measurements  of 
fundamental  mode  Rayleigh  wave  phase  velocity  have  extended  only  to  ~50  sec  period  (Hwang  and  Mitchell, 
1987;  Bhattacharya,  1992).  This  increased  period  range  of  our  dispersion  measurement  is  important  in 
constraining  the  deep  structure  of  the  south  Indian  Shield.  For  comparison,  the  dispersion  curve  for  the 
Canadian  Shield  (CANSD  -  Brune  and  Dorman,  1963)  is  also  plotted  in  Figure  8.  The  crustal  thickness  of  the 
south  Indian  Shield  from  the  receiver  function  study  discussed  below  is  36±1  km,  and  the  crust  of  model 
CANSD  is  35  km  thick.  The  short  period  (<35  sec)  and  long  period  (>130  sec)  portions  of  the  two  dispersion 
curves  are  similar,  but  between  40  and  125  sec  period  the  phase  velocities  for  the  south  Indian  Shield  are  low 
compared  to  those  of  model  CANSD.  The  preliminary  inversion  of  the  Indian  phase  velocity  curves  indicated 
the  presence  of  a  high  velocity  upper  mantle  lid  with  S  wave  velocity  ranging  from  -4.62  km/s  at  40-60  km 
depth  to  -4.72  km/s  at  80-120  km  depth  and  a  low  velocity  of  -4.5  km/s  at  150-220  km  depth. 

CONCLUSIONS  AND  RECOMMENDATIONS 


During  the  past  year,  we  have  initiated  a  study  to  collect  seismic  data  from  southern  Asia  for  the  purpose  of 
developing  high  quality  velocity  models  to  aid  in  regional  seismic  event  location.  We  have  installed  five 
broadband  sensors  in  southern  Asia  and  are  now  collecting  data.  Using  these  data  collected  from  the  current 
deployment,  along  with  other  data  sources  from  the  region,  we  have  completed  surface  wave  and  receiver 
function  studies.  Based  upon  surface  wave  phase  velocity  measurements,  we  have  developed  a  shear  wave 
velocity  structure  for  the  crust  and  upper  mantle  of  the  south  Indian  shield.  The  average  crustal 
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Figure  7.  Moho  depth  from  receiver  function  analysis. 


Period  (sec) 


Figure  8.  Long-period  fundamental  mode  Rayleigh  wave  phase  velocity  dispersion.  Left:  Match  of  the  far 

waveform  (solid  line)  and  near  seismogram  (dashed  line)  after  being  fdtered  by  the  determined  transfer 
function.  Right:  Upper  -  dispersion  curve  for  the  south  Indian  Shield  (dot-dash  line)  and  ±la  bounds 
(dash  lines)  and  dispersion  curve  for  model  CANSD  of  Brune  and  Dorman  (1963)  (light  solid  line); 
Lower  -  coherency  vs.  period. 
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structure  beneath  the  profile  consists  of  two  layers  (hi  12  km,  Vsl  3.65  km/s,  h2  23  km,  Vs2  3.81  km/s)  overlying 
a  mantle  with  S„  velocity  4.61  km/s.  Inversions  of  the  phase  velocity  curves  indicate  the  presence  of  a  high 
velocity  upper  mantle  lid  with  S  wave  velocity  ranging  from  -4.62  km/s  at  40-60  km  depth  to  -4.72  km/s  at  80- 
120  km  depth  and  a  low  velocity  of  -4.5  km/s  at  150-220  km  depth.  We  have  also  analyzed  receiver  function 
data  from  14  broadband  sites  in  southern  India  and  determine  the  velocity  structure  is  similar  to  the  surface 
wave  results.  The  Moho  is  at  almost  constant  depth  along  the  western  edge  of  the  eastern  Dharwar  craton. 
However,  the  two  sites  on  the  western  Dharwar  craton  show  the  Moho  there  is  about  5  km  deeper,  a  result  also 
noted  in  the  Deep  Seismic  Soundings  from  the  region.  To  the  south  of  the  shield  in  the  granulite  terrain  the  crust 
is  43  km  thick  We  have  also  modeled  the  regional  waveforms  of  a  moderate  earthquake  that  occurred  near 
Koyna,  India,  in  September  2000  and  of  the  26  January  2001  Bhuj  main  shock  and  its  aftershocks  to  calibrate 
paths  to  the  regional  stations.  In  the  next  stage  of  this  project,  we  will  continue  the  data  collection  and  velocity 
model  development,  however  our  research  focus  will  be  shifted  northward  to  central  and  eastern  India, 
including  the  Shillong  Plateau. 
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